home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zgemm.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  14.5 KB  |  460 lines
  1. *
  2. ************************************************************************
  3. *
  4. *     File of the COMPLEX*16       Level-3 BLAS.
  5. *     ==========================================
  6. *
  7. *     SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
  8. *    $                   BETA, C, LDC )
  9. *
  10. *     SUBROUTINE ZSYMM ( SIDE,   UPLO,   M, N,    ALPHA, A, LDA, B, LDB,
  11. *    $                   BETA, C, LDC )
  12. *
  13. *     SUBROUTINE ZHEMM ( SIDE,   UPLO,   M, N,    ALPHA, A, LDA, B, LDB,
  14. *    $                   BETA, C, LDC )
  15. *
  16. *     SUBROUTINE ZSYRK ( UPLO,   TRANS,     N, K, ALPHA, A, LDA,
  17. *    $                   BETA, C, LDC )
  18. *
  19. *     SUBROUTINE ZHERK ( UPLO,   TRANS,     N, K, ALPHA, A, LDA,
  20. *    $                   BETA, C, LDC )
  21. *
  22. *     SUBROUTINE ZSYR2K( UPLO,   TRANS,     N, K, ALPHA, A, LDA, B, LDB,
  23. *    $                   BETA, C, LDC )
  24. *
  25. *     SUBROUTINE ZHER2K( UPLO,   TRANS,     N, K, ALPHA, A, LDA, B, LDB,
  26. *    $                   BETA, C, LDC )
  27. *
  28. *     SUBROUTINE ZTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  29. *    $                   B, LDB )
  30. *
  31. *     SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  32. *    $                   B, LDB )
  33. *
  34. *     See:
  35. *
  36. *        Dongarra J. J.,   Du Croz J. J.,   Duff I.  and   Hammarling S.
  37. *        A set of  Level 3  Basic Linear Algebra Subprograms.  Technical
  38. *        Memorandum No.88 (Revision 1), Mathematics and Computer Science
  39. *        Division,  Argonne National Laboratory, 9700 South Cass Avenue,
  40. *        Argonne, Illinois 60439.
  41. *
  42. *
  43. ************************************************************************
  44. *
  45.       SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
  46.      $                   BETA, C, LDC )
  47. *     .. Scalar Arguments ..
  48.       CHARACTER*1        TRANSA, TRANSB
  49.       INTEGER            M, N, K, LDA, LDB, LDC
  50.       COMPLEX*16         ALPHA, BETA
  51. *     .. Array Arguments ..
  52.       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
  53. *     ..
  54. *
  55. *  Purpose
  56. *  =======
  57. *
  58. *  ZGEMM  performs one of the matrix-matrix operations
  59. *
  60. *     C := alpha*op( A )*op( B ) + beta*C,
  61. *
  62. *  where  op( X ) is one of
  63. *
  64. *     op( X ) = X   or   op( X ) = X'   or   op( X ) = conjg( X' ),
  65. *
  66. *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
  67. *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
  68. *
  69. *  Parameters
  70. *  ==========
  71. *
  72. *  TRANSA - CHARACTER*1.
  73. *           On entry, TRANSA specifies the form of op( A ) to be used in
  74. *           the matrix multiplication as follows:
  75. *
  76. *              TRANSA = 'N' or 'n',  op( A ) = A.
  77. *
  78. *              TRANSA = 'T' or 't',  op( A ) = A'.
  79. *
  80. *              TRANSA = 'C' or 'c',  op( A ) = conjg( A' ).
  81. *
  82. *           Unchanged on exit.
  83. *
  84. *  TRANSB - CHARACTER*1.
  85. *           On entry, TRANSB specifies the form of op( B ) to be used in
  86. *           the matrix multiplication as follows:
  87. *
  88. *              TRANSB = 'N' or 'n',  op( B ) = B.
  89. *
  90. *              TRANSB = 'T' or 't',  op( B ) = B'.
  91. *
  92. *              TRANSB = 'C' or 'c',  op( B ) = conjg( B' ).
  93. *
  94. *           Unchanged on exit.
  95. *
  96. *  M      - INTEGER.
  97. *           On entry,  M  specifies  the number  of rows  of the  matrix
  98. *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
  99. *           Unchanged on exit.
  100. *
  101. *  N      - INTEGER.
  102. *           On entry,  N  specifies the number  of columns of the matrix
  103. *           op( B ) and the number of columns of the matrix C. N must be
  104. *           at least zero.
  105. *           Unchanged on exit.
  106. *
  107. *  K      - INTEGER.
  108. *           On entry,  K  specifies  the number of columns of the matrix
  109. *           op( A ) and the number of rows of the matrix op( B ). K must
  110. *           be at least  zero.
  111. *           Unchanged on exit.
  112. *
  113. *  ALPHA  - COMPLEX*16      .
  114. *           On entry, ALPHA specifies the scalar alpha.
  115. *           Unchanged on exit.
  116. *
  117. *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  118. *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
  119. *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
  120. *           part of the array  A  must contain the matrix  A,  otherwise
  121. *           the leading  k by m  part of the array  A  must contain  the
  122. *           matrix A.
  123. *           Unchanged on exit.
  124. *
  125. *  LDA    - INTEGER.
  126. *           On entry, LDA specifies the first dimension of A as declared
  127. *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
  128. *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
  129. *           least  max( 1, k ).
  130. *           Unchanged on exit.
  131. *
  132. *  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
  133. *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
  134. *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
  135. *           part of the array  B  must contain the matrix  B,  otherwise
  136. *           the leading  n by k  part of the array  B  must contain  the
  137. *           matrix B.
  138. *           Unchanged on exit.
  139. *
  140. *  LDB    - INTEGER.
  141. *           On entry, LDB specifies the first dimension of B as declared
  142. *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
  143. *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
  144. *           least  max( 1, n ).
  145. *           Unchanged on exit.
  146. *
  147. *  BETA   - COMPLEX*16      .
  148. *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
  149. *           supplied as zero then C need not be set on input.
  150. *           Unchanged on exit.
  151. *
  152. *  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
  153. *           Before entry, the leading  m by n  part of the array  C must
  154. *           contain the matrix  C,  except when  beta  is zero, in which
  155. *           case C need not be set on entry.
  156. *           On exit, the array  C  is overwritten by the  m by n  matrix
  157. *           ( alpha*op( A )*op( B ) + beta*C ).
  158. *
  159. *  LDC    - INTEGER.
  160. *           On entry, LDC specifies the first dimension of C as declared
  161. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  162. *           max( 1, m ).
  163. *           Unchanged on exit.
  164. *
  165. *
  166. *  Level 3 Blas routine.
  167. *
  168. *  -- Written on 8-February-1989.
  169. *     Jack Dongarra, Argonne National Laboratory.
  170. *     Iain Duff, AERE Harwell.
  171. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  172. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  173. *
  174. *
  175. *     .. External Functions ..
  176.       LOGICAL            LSAME
  177.       EXTERNAL           LSAME
  178. *     .. External Subroutines ..
  179.       EXTERNAL           XERBLA
  180. *     .. Intrinsic Functions ..
  181.       INTRINSIC          DCONJG, MAX
  182. *     .. Local Scalars ..
  183.       LOGICAL            CONJA, CONJB, NOTA, NOTB
  184.       INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
  185.       COMPLEX*16         TEMP
  186. *     .. Parameters ..
  187.       COMPLEX*16         ONE
  188.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  189.       COMPLEX*16         ZERO
  190.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  191. *     ..
  192. *     .. Executable Statements ..
  193. *
  194. *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
  195. *     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
  196. *     B  respectively are to be  transposed but  not conjugated  and set
  197. *     NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A
  198. *     and the number of rows of  B  respectively.
  199. *
  200.       NOTA  = LSAME( TRANSA, 'N' )
  201.       NOTB  = LSAME( TRANSB, 'N' )
  202.       CONJA = LSAME( TRANSA, 'C' )
  203.       CONJB = LSAME( TRANSB, 'C' )
  204.       IF( NOTA )THEN
  205.          NROWA = M
  206.          NCOLA = K
  207.       ELSE
  208.          NROWA = K
  209.          NCOLA = M
  210.       END IF
  211.       IF( NOTB )THEN
  212.          NROWB = K
  213.       ELSE
  214.          NROWB = N
  215.       END IF
  216. *
  217. *     Test the input parameters.
  218. *
  219.       INFO = 0
  220.       IF(      ( .NOT.NOTA                 ).AND.
  221.      $         ( .NOT.CONJA                ).AND.
  222.      $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
  223.          INFO = 1
  224.       ELSE IF( ( .NOT.NOTB                 ).AND.
  225.      $         ( .NOT.CONJB                ).AND.
  226.      $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
  227.          INFO = 2
  228.       ELSE IF( M  .LT.0               )THEN
  229.          INFO = 3
  230.       ELSE IF( N  .LT.0               )THEN
  231.          INFO = 4
  232.       ELSE IF( K  .LT.0               )THEN
  233.          INFO = 5
  234.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  235.          INFO = 8
  236.       ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
  237.          INFO = 10
  238.       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
  239.          INFO = 13
  240.       END IF
  241.       IF( INFO.NE.0 )THEN
  242.          CALL XERBLA( 'ZGEMM ', INFO )
  243.          RETURN
  244.       END IF
  245. *
  246. *     Quick return if possible.
  247. *
  248.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  249.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  250.      $   RETURN
  251. *
  252. *     And when  alpha.eq.zero.
  253. *
  254.       IF( ALPHA.EQ.ZERO )THEN
  255.          IF( BETA.EQ.ZERO )THEN
  256.             DO 20, J = 1, N
  257.                DO 10, I = 1, M
  258.                   C( I, J ) = ZERO
  259.    10          CONTINUE
  260.    20       CONTINUE
  261.          ELSE
  262.             DO 40, J = 1, N
  263.                DO 30, I = 1, M
  264.                   C( I, J ) = BETA*C( I, J )
  265.    30          CONTINUE
  266.    40       CONTINUE
  267.          END IF
  268.          RETURN
  269.       END IF
  270. *
  271. *     Start the operations.
  272. *
  273.       IF( NOTB )THEN
  274.          IF( NOTA )THEN
  275. *
  276. *           Form  C := alpha*A*B + beta*C.
  277. *
  278.             DO 90, J = 1, N
  279.                IF( BETA.EQ.ZERO )THEN
  280.                   DO 50, I = 1, M
  281.                      C( I, J ) = ZERO
  282.    50             CONTINUE
  283.                ELSE IF( BETA.NE.ONE )THEN
  284.                   DO 60, I = 1, M
  285.                      C( I, J ) = BETA*C( I, J )
  286.    60             CONTINUE
  287.                END IF
  288.                DO 80, L = 1, K
  289.                   IF( B( L, J ).NE.ZERO )THEN
  290.                      TEMP = ALPHA*B( L, J )
  291.                      DO 70, I = 1, M
  292.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  293.    70                CONTINUE
  294.                   END IF
  295.    80          CONTINUE
  296.    90       CONTINUE
  297.          ELSE IF( CONJA )THEN
  298. *
  299. *           Form  C := alpha*conjg( A' )*B + beta*C.
  300. *
  301.             DO 120, J = 1, N
  302.                DO 110, I = 1, M
  303.                   TEMP = ZERO
  304.                   DO 100, L = 1, K
  305.                      TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J )
  306.   100             CONTINUE
  307.                   IF( BETA.EQ.ZERO )THEN
  308.                      C( I, J ) = ALPHA*TEMP
  309.                   ELSE
  310.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  311.                   END IF
  312.   110          CONTINUE
  313.   120       CONTINUE
  314.          ELSE
  315. *
  316. *           Form  C := alpha*A'*B + beta*C
  317. *
  318.             DO 150, J = 1, N
  319.                DO 140, I = 1, M
  320.                   TEMP = ZERO
  321.                   DO 130, L = 1, K
  322.                      TEMP = TEMP + A( L, I )*B( L, J )
  323.   130             CONTINUE
  324.                   IF( BETA.EQ.ZERO )THEN
  325.                      C( I, J ) = ALPHA*TEMP
  326.                   ELSE
  327.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  328.                   END IF
  329.   140          CONTINUE
  330.   150       CONTINUE
  331.          END IF
  332.       ELSE IF( NOTA )THEN
  333.          IF( CONJB )THEN
  334. *
  335. *           Form  C := alpha*A*conjg( B' ) + beta*C.
  336. *
  337.             DO 200, J = 1, N
  338.                IF( BETA.EQ.ZERO )THEN
  339.                   DO 160, I = 1, M
  340.                      C( I, J ) = ZERO
  341.   160             CONTINUE
  342.                ELSE IF( BETA.NE.ONE )THEN
  343.                   DO 170, I = 1, M
  344.                      C( I, J ) = BETA*C( I, J )
  345.   170             CONTINUE
  346.                END IF
  347.                DO 190, L = 1, K
  348.                   IF( B( J, L ).NE.ZERO )THEN
  349.                      TEMP = ALPHA*DCONJG( B( J, L ) )
  350.                      DO 180, I = 1, M
  351.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  352.   180                CONTINUE
  353.                   END IF
  354.   190          CONTINUE
  355.   200       CONTINUE
  356.          ELSE
  357. *
  358. *           Form  C := alpha*A*B'          + beta*C
  359. *
  360.             DO 250, J = 1, N
  361.                IF( BETA.EQ.ZERO )THEN
  362.                   DO 210, I = 1, M
  363.                      C( I, J ) = ZERO
  364.   210             CONTINUE
  365.                ELSE IF( BETA.NE.ONE )THEN
  366.                   DO 220, I = 1, M
  367.                      C( I, J ) = BETA*C( I, J )
  368.   220             CONTINUE
  369.                END IF
  370.                DO 240, L = 1, K
  371.                   IF( B( J, L ).NE.ZERO )THEN
  372.                      TEMP = ALPHA*B( J, L )
  373.                      DO 230, I = 1, M
  374.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  375.   230                CONTINUE
  376.                   END IF
  377.   240          CONTINUE
  378.   250       CONTINUE
  379.          END IF
  380.       ELSE IF( CONJA )THEN
  381.          IF( CONJB )THEN
  382. *
  383. *           Form  C := alpha*conjg( A' )*conjg( B' ) + beta*C.
  384. *
  385.             DO 280, J = 1, N
  386.                DO 270, I = 1, M
  387.                   TEMP = ZERO
  388.                   DO 260, L = 1, K
  389.                      TEMP = TEMP +
  390.      $                      DCONJG( A( L, I ) )*DCONJG( B( J, L ) )
  391.   260             CONTINUE
  392.                   IF( BETA.EQ.ZERO )THEN
  393.                      C( I, J ) = ALPHA*TEMP
  394.                   ELSE
  395.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  396.                   END IF
  397.   270          CONTINUE
  398.   280       CONTINUE
  399.          ELSE
  400. *
  401. *           Form  C := alpha*conjg( A' )*B' + beta*C
  402. *
  403.             DO 310, J = 1, N
  404.                DO 300, I = 1, M
  405.                   TEMP = ZERO
  406.                   DO 290, L = 1, K
  407.                      TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L )
  408.   290             CONTINUE
  409.                   IF( BETA.EQ.ZERO )THEN
  410.                      C( I, J ) = ALPHA*TEMP
  411.                   ELSE
  412.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  413.                   END IF
  414.   300          CONTINUE
  415.   310       CONTINUE
  416.          END IF
  417.       ELSE
  418.          IF( CONJB )THEN
  419. *
  420. *           Form  C := alpha*A'*conjg( B' ) + beta*C
  421. *
  422.             DO 340, J = 1, N
  423.                DO 330, I = 1, M
  424.                   TEMP = ZERO
  425.                   DO 320, L = 1, K
  426.                      TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) )
  427.   320             CONTINUE
  428.                   IF( BETA.EQ.ZERO )THEN
  429.                      C( I, J ) = ALPHA*TEMP
  430.                   ELSE
  431.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  432.                   END IF
  433.   330          CONTINUE
  434.   340       CONTINUE
  435.          ELSE
  436. *
  437. *           Form  C := alpha*A'*B' + beta*C
  438. *
  439.             DO 370, J = 1, N
  440.                DO 360, I = 1, M
  441.                   TEMP = ZERO
  442.                   DO 350, L = 1, K
  443.                      TEMP = TEMP + A( L, I )*B( J, L )
  444.   350             CONTINUE
  445.                   IF( BETA.EQ.ZERO )THEN
  446.                      C( I, J ) = ALPHA*TEMP
  447.                   ELSE
  448.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  449.                   END IF
  450.   360          CONTINUE
  451.   370       CONTINUE
  452.          END IF
  453.       END IF
  454. *
  455.       RETURN
  456. *
  457. *     End of ZGEMM .
  458. *
  459.       END
  460.